clear; clc; close all;

k = 3;
h = 0.01;
T = 600;
transient = 550;
N = round(T/h);
N_trans = round(transient/h);

mu_vals = linspace(-3*pi, 3*pi, 100);
eta_vals = linspace(-3*pi, 3*pi, 100);

SE_map = zeros(length(eta_vals), length(mu_vals));

parfor i_eta = 1:length(eta_vals)
    eta = eta_vals(i_eta);
    local_SE_row = zeros(1, length(mu_vals));
    for i_mu = 1:length(mu_vals)
        mu = mu_vals(i_mu);
        
        x = [1e-9; 0; 0; mu; eta];
        f = @(x) system_ode(x, k);

        % 丢弃暂态
        for t = 1:N_trans
            x = RK4_step(f, x, h);
        end

        % 采集轨迹y变量
        y_seq = zeros(N - N_trans, 1);
        for t = 1:(N - N_trans)
            x = RK4_step(f, x, h);
            y_seq(t) = x(2); % y变量
        end

        % 计算谱熵
        local_SE_row(i_mu) = spectral_entropy(y_seq);
    end
    SE_map(i_eta, :) = local_SE_row;
end

% 绘图
figure;
imagesc(mu_vals/pi, eta_vals/pi, SE_map);
set(gca,'YDir','normal');

colormap(jet);
colorbar;
caxis([0 0.5]); % 谱熵归一化范围0-1

xlabel('\mu');
ylabel('\eta');

xticks(-3:1:3);
yticks(-3:1:3);

xticklabels({'-3\pi','-2\pi','-\pi','0','\pi','2\pi','3\pi'});
yticklabels({'-3\pi','-2\pi','-\pi','0','\pi','2\pi','3\pi'});

grid on;

%% 子函数：系统方程
function dx = system_ode(x, k)
    dx = zeros(5,1);
    dx(1) = x(2) + x(3) - k * cos(x(5)) * x(2);
    dx(2) = -x(1) + x(3);
    dx(3) = -x(1) - x(3) + k * cos(x(4)) * x(1);
    dx(4) = x(1);
    dx(5) = x(2);
end

%% RK4一步积分
function x_next = RK4_step(f, x, h)
    k1 = f(x);
    k2 = f(x + h/2 * k1);
    k3 = f(x + h/2 * k2);
    k4 = f(x + h * k3);
    x_next = x + h/6*(k1 + 2*k2 + 2*k3 + k4);
end

%% 谱熵计算函数
function SE = spectral_entropy(x)
    N = length(x);
    Pxx = abs(fft(x)).^2;
    Pxx = Pxx(1:floor(N/2));
    Pxx = Pxx / sum(Pxx);
    Pxx(Pxx == 0) = [];
    SE = -sum(Pxx .* log2(Pxx)) / log2(length(Pxx));
end
